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ABSTRACT 

We use global magnetohydrodynamic simulations, in a pseudo-Newtonian potential, 
to investigate the temporal variability of accretion discs around Schwarzschild black 
holes. We use the vertically-averaged magnetic stress in the simulated disc as a proxy 
for the rest-frame dissipation, and compute the observed emission by folding this 
through the transfer function describing the relativistic beaming, light bending and 
time delays near a non-rotating black hole. The temporal power spectrum of the pre- 
dicted emission from individual annuli in the disc is described by a broken power law, 
with indices of « —3.5 at high frequency and « to —1 at low frequency. Integrated 
over the disc, the power spectrum is approximated by a single power law with an 
index of —2. Increasing inclination boosts the relative power at frequencies around 
~ 0.3/„is, where /ms is the orbital frequency at the marginally stable orbit, but no 
evidence is found for sharp quasi-periodic oscillations in the lightcurve. Assuming that 
fluorescent iron line emission locally tracks the continuum flux, we compute simulated 
broad iron line profiles. We find that relativistic beaming of the non-axisymmetric 
emission profile, induced by turbulence, produces high-amplitude variability in the 
iron line profile. We show that this substructure within the broad iron line profile 
can survive averaging over a number of orbital periods, and discuss the origin of the 
anomalous X-ray spectral features, recently reported by Turner et al. (2002) for the 
Seyfert galaxy NGC 3516, in the context of turbulent disc models. 

Key vifords: accretion, accretion discs — black hole physics — MHD — turbulence 
— X-rays: binaries — galaxies: active 



1 INTRODUCTION 

A major aim of the astrophysical study of black holes is to 
determine the fundamental parameters of a rotating black 
hole - the mass M and spin parameter a - from astronomi- 
cal observations. For Galactic black hole candidates such as 
Cygnus X-1, GRO J1655-40 and GRS 1915+105, the tem- 
poral power spectrum holds promise as a route to achiev- 
ing this goal (van der Klis 2000). The variability of these 
sources is conventionally described as the superposition of a 
broadband noise component with one or more quasi-periodic 
oscillations (QPOs), whose frequencies presumably encode 
information about the inner accretion flow and black hole 
properties. For AGN, where the long dynamical times asso- 
ciated with supermassive black holes hamper similar stud- 
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ies, the main focus is on inferring black hole properties from 
the profiles of relativistically broadened iron emission lines 
arising from near the marginally stable orbit (Fabian et al. 
1989; Tanaka et al. 1995; Nandra et al. 1997; Wilms et al. 
2001; Lee et al. 2002). Iron line profiles themselves are highly 
variable (Iwasawa et al. 1996). 

For accreting black holes, and for other relativistically 
compact objects, the observed lightcurve arises from the 
modulation of rest-frame variability (i.e. variability that 
would be detected by an observer orbiting with the disc 
gas) by photon propagation elTects. Since there is now a 
consensus that magnetorotational instabilities (MRI, Balbus 
& Hawley 1991; Balbus & Hawley 1998) provide a source 
of angular momentum transport in well ionized discs, an 
unavoidable source of rest-frame variability is the fluctuat- 
ing dissipation inherent to magnetohydrodynamic (MHD) 
disc turbulence. Simulations of the MRI show that the high 
frequency temporal power spectra of the mass accretion 
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rate and magnetic stress display steep power-law spectra 
(Kawaguchi et al. 2000; Hawley & Krolik 2001, 2002), and 
suggest that the broadband noise component of observed 
power density spectra may be identified with the temporal 
fluctuations of MHD disc turbulence. Spatial fluctuations in 
the dissipation may also play an important role. For a disc 
inclined to the line of sight, beaming enhances the emission 
from the approaching side of the disc relative to that from 
the receeding side. The combination of this photon propaga- 
tion effect with any source of persistent non-axisymmetric 
dissipation in the inner regions of the disc could lead to 
high frequency QPOs in the lightcurve (e.g. Abramowicz et 
al. 1991; Karas 1999). 

In this paper, we use non-rclativistic global MHD sim- 
ulations of the MRI (Armitage 1998; Hawley 2000; Arlt & 
Rudiger 2001) to model the rest-frame variability of accre- 
tion onto a Schwarzschild black hole. A pseudo-Newtonian 
potential (Paczynski & Wiita 1980) is used to describe the 
gravitational potential in the vicinity of such a black hole. 
We extend recent three-dimensional numerical studies of 
black hole accretion (Hawley 2000; Hawley & Krolik 2001; 
Armitage, Reynolds & Chiang 2001; Hawley 2001; Reynolds 
& Armitage 2001; Hawley, Balbus & Stone 2001; Machida, 
Matsumoto & Mineshige 2001; ; Hawley & Krohk 2002; Igu- 
menshchev & Narayan 2002; Machida & Matsumoto 2003; 
Igumenshchev, Narayan & Abramowicz 2003) by taking full 
account of the distinction between the rest-frame disc emis- 
sion and that seen by a distant observer. We use these sim- 
ulations to study both the predicted power spectra from 
black holes accreting via a geometrically thin disc, and the 
expected variability of the iron line emission that has been 
observed in the X-ray spectra of Seyfcrt galaxies (Tanaka 
et al. 1995; Wilms et al 2001; Lee et al. 2002) and Galactic 
Black Hole Candidates (Balucinska-Church & Church 2000; 
Miller et al. 2002). The details of the MHD simulations are 
presented in Section 2, while Section 3 describes our com- 
putation of the observed continuum and line emission tak- 
ing into account relativistic photon propagation. Predicted 
temporal power spectra are presented in Section 4. In Sec- 
tion 5, we explore our predictions for variability of the broad 
iron line. It is suggested that the non-axisymmetric structure 
discussed in this paper may already have been observed in 
Chandra High Energy Transmission Grating (HETG) data 
for the Seyfert galaxy NGC 3516, and we proceed to quan- 
tify the nature of the variability that might be observed 
in the future, using higher signal-to-noiso Constellation-X 
data. Our results, predictions and speculations are summa- 
rized in Section 6. 



2 THE MHD SIMULATIONS 

We model the accretion disc around a non-rotating black 

hole using ideal, Newtonian MHD. The fluid is assumed to 
be isothermal, with sound speed Cs, so the equations to be 
solved are. 



(4) 
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where S = pv is the momentum field and the other symbols 
have their usual meanings. These equations are solved us- 
ing the ZEUS MHD code (Stone & Norman 1992a, 1992b), 
operating in cylindrical polar co-ordinates {z,r,4>). ZEUS is 
a fixed-grid, time-explicit Eulerian code which uses an arti- 
ficial viscosity to handle shocks, and a combination of con- 
strained transport (Evans & Hawley 1988) and the method 
of characteristics (Stone & Norman 1992b) to evolve the 
magnetic fields. When operated using van Leer advection, 
as in this work, it is formally of second order accuracy. 

A pseudo-Newtonian gravitational potential (Paczynski 
& Wiita 1980) is used to mimic the effects of general rela- 
tivity near a non-rotating black hole of mass M, 

$ = , 5) 

r-rg 

where r-g = 2GM/c^. This potential reproduces the pres- 
ence of a marginally stable circular orbit at rms = 6GM/c^ . 
For computational reasons, we also omit the vertical compo- 
nent of the gravitational force. This 'cylindrical disc' approx- 
imation has been used in several previous studies (Hawley, 
Gammie & Balbus 1995; Armitage 1998; Reynolds & Ar- 
mitage 2001; Hawley 2001; Steinacker & Papaloizou 2002), 
and greatly reduces the computational cost of global disc 
simulations. Of course, it also precludes study of physical 
effects such as the Parker instability and the development 
of a disc corona (e.g. Miller & Stone 2000). 

To investigate variability caused by the combination of 
a turbulent disc with relativistic beaming effects, it is nec- 
essary to model all 27r of the disc in azimuth. Apart from 
this (trivial) extension, our numerical setup is the same as 
that which we have used previously (Armitage, Reynolds & 
Chiang 2001; Reynolds & Armitage 2001). A hydrodynami- 
cally stable disc flow, threaded by a weak vertical magnetic 
field, is set up as the initial state of the simulations. We be- 
gin by defining a Gaussian surface density profile, centred at 
4 rms, where rms is the radius of the marginally stable circu- 
lar orbit, with a width of 2.66 rms. This disc, initialized with 
Keplerian angular velocity and zero magnetic fields, is run 
with the code in ID mode until a numerical equilibrium is 
attained. We then add a weak vertical magnetic field (ratio 
of gas to magnetic pressure /3 — 10'') at all radii where the 
density exceeds a threshold value (set to be approximately 
10 percent of the max;imum density), and continue the run 
in 3D. The boundary conditions are periodic in both z and 
0, and set to outflow at the radial boundaries. 'Outflow' 
boundary conditions are implemented in ZEUS by setting 
the values of all physical quantities in boundary zones to be 
equal to their values in the active zones at the edge of the 
grid. 

Two simulations were run, differing only in the adopted 
sound speed. One run had a ratio of the sound speed to 
the Keplerian velocity at rms of Cs/v^ = 0.061, while a sec- 
ond run had Cs/v^ = 0.041 at rms- We will refer to these 
as the 'hot' and 'cold' run respectively, though the actual 
temperature differences are modest. Both simulations were 
evolved for a total of 140 orbital periods at rms, in identical 
computational domains. 
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Figure 1. The disc mass as a function of time for tlie 'liot' (solid 
curve and 'cold' (dashed curve) runs, shown as a fraction of the 
initial disc mass. Only the fully turbulent portion of the run with 
t > 50, where t is measured in units of the orbital period at rms, 
is used for analysis of the variability. 

< <j) <2tt, (8) 

with 48 grid points in z, 192 grid points in r, and 288 grid 
points in </!>. A non-uniform radial grid, with mesh points 
spaced such that r,;+i = (1 + S)ri, with S a constant, is 
employed, so that Ar/A<^ is constant at all radii. The inner 
radial boundary was placed as close to the marginally stable 
orbit as possible (to maximise the timestep which is almost 
always determined by the azimuthal velocity at rin), but far 
enough within Vms that the flow across the boundary was 
both supersonic and super- Alfvenic. As has been noted by 
Igumenshchev et al. (2003), this inner boundary condition 
becomes inappropriate if large net amounts of magnetic flux 
are able to be dragged inward by the accretion flow. Here, 
we make the assumption - justified in part by the resuhs of 
Lubow, Papaloizou & Pringle (1994) - that extensive field 
dragging does not occur for thin discs. 

2.1 Properties of the simulations 

The basic properties of the MHD disc turbulence seen in 
the simulations are, unsurprisingly, very similar to those re- 
ported previously (i.e. although we need a full 2tt in azimuth 
to study variability, smaller azimuthal domains suffice for 
most other purposes). Here, we summarize a few features 
that are pertinent to the subsequent analysis. 

Figure0shows the total mass in the computational do- 
main as a function of time for the hot and cool runs re- 
spectively. Significant accretion commences after around 30 
orbits of evolution at the marginally stable orbit. For the hot 
run, the accretion rate is very roughly constant after around 
50 orbits of evolution, and it is this portion of the simulation 
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Figure 2. The vertically integrated magnetic stress for the hot 
disc at a radius of 1.5 rms, evaluated at a single azimuthal point 
corotating with the local Keplerian angular velocity. The upper 
panel shows the raw measurement, the lower panel the same curve 
smoothed over an interval equal to the local orbital period. The 
stress is shown normalized to the instantaneous average value at 
the same radius. 

which we analyse for variability. Note, however, that there 
remains a secular trend in the accretion rate even over this 
latter portion (initially increasing, then decreasing once a 
significant fraction of the initial mass has been accreted). 
Since this is an artefact of the computational setup, we can 
only look at variability on substantially shorter timescales 
than the run length^ 

A single persistent hotspot, orbiting near the innermost 
stable orbit where relativistic effects are most pronounced, 
would lead to a strong periodic feature in the lightcurve 
(Karas 1999). Such QPOs will be diluted in more realistic 
models, first by the existence of patchy dissipation at all 
radii, and second by the limited lifetime of regions of en- 
hanced dissipation. Figure |21 illustrates the latter point. We 
plot the magnetic stress, vertically integrated through the 
disc, at a single point corotating with the local Keplerian 
orbital velocity. Results for the hot simulation are shown, 
although very similar results are obtained for the cold case. 
Large fluctuations are observed, with the peak stress exceed- 
ing the local mean by a factor that can be as large as ~ 4. 
However, these are also transient features, with a tempo- 
ral coherence time of only a few orbital periods. If, as we 
assume later, the magnetic stress provides a proxy for the 

^ Because we have only one realization of each simulation, there 
are also large uncertainties in our estimation of the temporal 
power spectrum at the lowest frequencies. These purely statistical 
uncertainties would preclude us from making statements about 
variability on timescales comparable to the run length, even if 
the mean accretion rate showed no secular trend. 
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Figure 3. The magnetic torque at the marginally stable orbit as 
a fraction of the torque in the disc at a radius of 2rms for the 
simulations with high (solid line) and low (dashed line) sound 
speed. Non-zero stress at rms becomes relatively more important 
for the disc structure as the disc temperature increases. 



disc dissipation, it is clear that 'hotspots' arising from MHD 
disc turbulence will lead to at most broad QPOs in power 
spectra of the disc emission. 

Recent analytic (Krolik 1999; Gammie 1999; Agol & 
Krolik 2000; Afshordi & Paczynski 2003) and numerical 
work (Hawley 2000; Armitage, Reynolds & Chiang 2001; 
Reynolds & Armitage 2001; Hawley & KroUk 2001, 2002) 
has investigated the possibility that the structure of the in- 
ner accretion disc might be modified by the existence of 
magnetic coupling to gas in the plunging region at r < rms- 
The presence of stress at the last stable orbit is potentially 
important for studies of variability, as a disc with non-zero 
stress at rms has a dissipation profile that is more centrally 
concentrated than a standard disc (Agol & Krolik 2000). 
Figure 13 shows how important this effect is for our simula- 
tions. We plot the mean stress at rms as a fraction of the 
stress at 2rms (this ratio is adopted because the changing 
accretion rate means that the stress at rms does not tend 
to a constant value). At late times, we find that this ra- 
tio fiuctuates around a value of around 0.3 for the hot disc 
run. There are, however, large fiuctuations, which persist 
over timescales that can be as long as 10 orbits at rms- As 
discussed by Papaloizou & Nelson (2003), and by Winters, 
Balbus & Hawley (2003), these persistent fluctuations are a 
generic feature of disc MHD turbulence. They are present 
in local (e.g. Brandenburg et al. 1995; Winters, Balbus & 
Hawley 2003) as well as global simulations (Papaloizou & 
Nelson 2003) , and their existence means that averaging over 
long timescales is needed to recover the mean properties of 
the flow. 

As in previous simulations (Reynolds & Armitage 
2001), there is evidence that the contribution of stress at 



rms becomes relatively less important as the sound speed 
(or, equivalently, the implied disc thickness h/r) is reduced. 
In our subsequent analysis, we compute both raw lightcurves 
(which include disc dissipation arising from stress at rms), 
and normalized lightcurves which explicitly fix the form 
of the radial dissipation profile to match the predictions 
of zero-torque boundary condition disc models. Comparing 
these lightcurves allows us to investigate whether the pres- 
ence of stress at the marginally stable orbit has a significant 
(and potentially observable) effect on the predicted tempo- 
ral power spectrum. 



3 RELATIVISTIC TRANSFER FUNCTION 

With a number of approximations, described in Section 4.1, 
the MHD simulations can be used to define a plausible 
time-dependent pattern of emission across the disc plane, 
F{r,(j),t). When the region of interest is close to the black 
hole, as it is in this study, one must also account for general 
relativistic effects (i.e., beaming, light bending and time- 
delays) in order to determine the appearance of the simu- 
lated disc to an observer at infinity. In this section, we briefiy 
discuss how we incorporate these effects. 

Suppose the observer is viewing the disc from a large 
distance at an inclination of i from the normal to the disc 
plane. We define cartesian coordinates {x,y) on the ob- 
server's image plane such that the origin x — y = lies 
at the center of the event horizon's image and the y-axis 
is parallel with the apparent (i.e. projected) normal to the 
disc. Let us consider a photon that is emitted from the disc 
at position (r, 0) with frequency and is subsequently ob- 
served, a Schwarzschild time t later, at position {x,y) on the 
image plane with frequency Vohs- Our task is to compute the 
functions x{r,4>-i), y{r,4>;i), t{x,y;i), and g{x,y;i), where 
we have defined the frequency shift factor as 

5(^,y) = — ■ (9) 

For a number of inclinations (i — 1°, 10°, 20°, 80°), 
we integrated null geodesies (photon paths) through the 
Schwarzschild metric 



ds 



1 



2GM 



dr 



+ r{de^+ sin 



(10) 



which exactly describes the space-time geometry around a 
non-rotating and uncharged black hole. Starting from a par- 
ticular point {x,y) on the image plane, we integrate photon 
paths down to the disc plane 9 — 71/2 using the first integrals 
of the geodesic equation (conservation of angular momentum 
and energy, arising from the (p and t independence of the 
metric, respectively), the Carter constant, and the light-like 
nature of the path. For these details, we refer the reader to 
Appendix A of Reynolds et al. (1999). Turning points in the 
photon orbits are identified and treated using the method of 
Ranch & Blandford (1994). These integrations allow us to 
determine the maps x{r,(f)\i), y{r,(f);i) and t{x,y;i). 

The frequency shift factor for a given photon path, g, 
is given by 



9 - 



(11) 
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where is the 4-momentum of the photon at the 
disc/observer, and u'^ is the 4- velocity of the disc/observer. 
The photon 4-momentum is determined by the photon inte- 
gration, and we suppose that the observer is stationary (in 
Schwarzschild coordinates) and very distant from the black 
hole; thus u'^^^^ = (1,0,0,0). Hence, we can compute g for 
any given photon path provided we know the 4- velocity field 
on the disc plane. Initially, one might be tempted to take the 
3-velocity field directly from our MIfD simulations (which 
are approximately 'Keplerian' with respect to the modified 
potential) and construct the appropriate 4-velocity field as 
input for our computation of g. However, this approach is 
problematic for the following reason. Since our MHD simu- 
lation is intrinsically non-relativistic, and based upon a po- 
tential that diverges at the event horizon, the 3-velocities 
within the simulated plunging region can significantly ex- 
ceed those that would be found in a relativistic calculation. 
This would lead to artificially high values of \g — 1|, and 
hence artificially strong relativistic beaming. 

Motivated by the fact that the velocity field in our sim- 
ulated disc is very close to Keplerian with respect to the 
pseudo-Newtonian potential, we choose instead to use the 
correct relativistic 4-velocity for test-particle circular mo- 
tion in our computation of g. This is well defined down to 
the photon circular orbit at 3GM/ c? , which lies inside the in- 
ner edge of our computational grid. Of course, disc material 
within some radius close to the radius of marginal stability 
(at QGM/c^) will no Ion ger be on circular orbits, but the 
deviations in the relevant parts of the disc are small. Em- 
ploying this assumption, we compute the function g{x,y; i). 



4 TEMPORAL VARIABILITY OF THE 
SIMULATED FLUX 

4.1 The simulated lightcurve 

A direct calculation of the emission from the disc would re- 
quire, first, that we solve an energy equation, and second 
that we attain sufficient resolution in a global simulation to 
model a stratified disc and corona (Miller & Stone 2000) 
where a significant fraction of the dissipation probably oc- 
curs. Since our actual calculation is more limited than this, 
we must instead choose some proxy for the dissipation. Out 
of several possible proxies, we have opted to follow Hawley & 
Krolik (2001), and use the magnetic stress as a local tracer 
of where in the fiow dissipation is expected to occur ^. We 
then compute lightcurves and power spectra by making the 
further assumption that the dissipated energy is promptly 
radiated from the disc. 

The most important assumption is obviously the afore- 
mentioned one - that at a given radius, the azimuthal de- 
pendence of the emission is proportional to the vertically in- 
tegrated magnetic stress. Even given this assumption, how- 
ever, there is some ambiguity as to the best way to calculate 
the lightcurve. As described below, we calculate both 'raw' 



^ Data storage limitations prompted us to make this basic choice 
prior to running full resolution simulations. During the runs, ver- 
tically integrated slices of BrB^jy were the only quantities output 
at high time resolution. 



lightcurves, and lightcurves in which the average disc dissi- 
pation has been normalized to match an analytic form. 



4-1.1 The raw lightcurve 

To compute the raw lightcurve, we first calculate the emis- 
sion map -Fraw that would be seen by observers orbiting with 
the disc fiow. fi-aw is given by (e.g. Hubeny & Hubeny 1998). 



-Fraw(r, <;/>,*) 



GM ^A\ 
r'^ \b) 



Br B^dz 



(12) 



where the relativistic correction factors are. 



A 
B 



= 1 



= 1 



2GM 



3GM 



(13) 



Using the transfer function, which gives the mapping be- 
tween (r, <j)) in the disc plane and {x, y) in the image plane, 
we then compute the image that would be seen by a dis- 
tant observer viewing the disc at a particular inclination 
angle, i. We assume that the observer images the disc in a 
fixed spectral band, which introduces a 'K-correction' that 
depends upon the spectrum of the disc emission. If at each 
radius the energy spectrum is oc (fairly typically for 
accreting black holes), the resulting image is. 



I{x, y, t) oc g^{x, y)Fr^„{r, 0, tcmit). 



(14) 



where one factor of g comes from the K-correction. We al- 
low for the varying flight time of photons from different 
parts of the accretion disc by distinguishing between the 
time temit, at which photons are emitted from the disc, and 
time t when the photons are observed. Examples of these 
images are shown in Figure|l| The raw lightcurve at time ti, 
J-{ti) is then simply the sum over the image of the intensity 
at time ti at each pixel {xj,yk), 



^{ti) = y^J{xj,yk,ti 



(15) 



We emphasize that the predicted variability properties de- 
pend upon the assumed emission spectrum from the disc as 
well as on the statistical fluctuations of MHD disc turbu- 
lence. In particular, modeling blackbody emission from the 
disc would require a different radial weighting of the simu- 
lated disc dissipation. 

The resulting light curves are shown in Figure |K| for 
several different inclinations. As already noted, the shape 
of these lightcurves over the longest timescales reflects only 
the viscous evolution of the ring of gas used for our initial 
conditions. There is a secular increase in the inner accretion 
rate (and accompanying dissipation) as the MRI saturates 
at increasingly large radii, followed (for the hot simulation) 
by an eventual decline once a substantial fraction of the gas 
has been accreted. These slow changes in accretion rate are 
artefacts of the initial conditions used in the simulation — 
real accretion discs are long lived objects with very large 
reservoirs of material. Of more relevance to real systems is 
the inclination dependence of these light curves. One ob- 
serves a marked increase in rapid variability for the more 
edge on systems. This is readily interpreted as the effects 
of relativistic beaming coupled with the 'blobby' nature of 
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Figure 4. View of the disc as seen by a distant observer at an inclination angle of 5° (upper left), 30° (upper right), 55° (lower 
left) and 80° (lower right). In these raw images, note the presence of stress extending to the inner boundary of the computa- 
tional domain, within the marginally stable circular orbit. Movies showing the evolution of the simulated disc are available at 
[http : // j llawwM . Col orado . edu/~pj a/black_hole .html| 



the simulated emission. We shall explore this point in more 
detail below when we discuss temporal power spectra. 



J^.1.2 Novikov-Thome normalized lightcurve 

The raw lightcurve computed from the simulation self- 
consistently includes the contribution from torque at the 
marginally stable orbit to the emission of the inner disc. As 
shown in Figure^ this contribution can be important, espe- 
cially for hotter (and therefore geometrically thicker) discs 
(Reynolds & Armitage 2001). However, because the domain 
of the simulation is limited, only the inner region of the sim- 
ulated disc is in a steady-state. The raw lightcurve therefore 
underestimates the contribution to the flux from larger radii 
in the disc. 

Cognisant of this drawback, we also compute lightcurves 
by combining the simulated fluctuations in the stress (in (j) 
and t) with an analytic dissipation profile due to Novikov & 
Thorne (1973). Specializing to the case of a Schwarzschild 
black hole, the radial dependence of the disc flux measured 
by an observer orbiting with the disc gas is (Page & Thorne 



1974), 

FNT(r) oc 

where. 



x'^{x''^ — 3x) 



r- \f?, ( X - \f?, 

X ~ Vb — In 



2 V^-^ 
, Vs. f x + V3_ 



{GM/c^ 



(16) 



(17) 



Given this dissipation profile, we proceed to construct a 
lightcurve as before, except that we replace the emission 
map Fraw (r, 0, t) with F]<iT{r, 4>,t), defined as, 

BrB4,{r,(l>,t) 



-Fnt (r, 0, t) = 



< BrB^ > (r) 



FNT(r) 



(18) 



In this expression Bj.B^ in both the numerator and the de- 
nominator is understood to be vertically integrated, and the 
angle brackets denote averaging over both </> and t. By con- 
struction, the resulting dissipation map F^T{r,(f),t) has the 
same time-averaged dissipation profile as a Novikov- Thorne 
disc, on top of which are superimposed fluctuations whose 
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Figure 5. Raw lightcurves showing the integrated flux from the disc as a function of the system inclination. From top downward, 
j = 1°, 20°, 40°, 60°, 80°. The small scale fluctuations increase with increasing system inclination. The hot simulation is displayed in the 
left panel, whereas the right panel displays the cold simulation. 



statistical properties are computed directly from the MHD 
simulation. Using this normalized dissipation map (which 
we shall refer to as a 'NT-normalization'), one can construct 
the observed emission map and a lightcurve in exactly the 
manner discussed above. 



4.2 The temporal power spectrum 

For a given lightcurve, we compute the temporal power spec- 
trum G{f), where / is the frequency, via a straightforward 
Fourier transform together with the "Leahy normalization" , 



G{.f) = 



(19) 



We truncate the resulting power spectra at the high- 
frequency end at the Nyquist frequency (i.e., half of the 
sampling frequency of the lightcurve) to account for the fi- 
nite sampling time. In addition, frequencies are rescaled such 
that the orbital frequency at the radius of marginal stability 
is unity. 

Figure |S| shows power spectra computed using emission 
that originates from several narrow annuli in the inner disc 
of the hot simulation. The power spectra are described well 
by broken power laws. Below a break frequency, which is 
comparable to the local orbital frequency, G{f) is flat or 
modestly declining {G{f) oc f~^). At higher frequencies, 
above the break, we obtain a steep spectrum with G{f) oc 

^-3.5 

Power spectra computed using the integrated disc emis- 
sion (ie directly from the lightcurves in Figure]^ are shown 
in Figures Q and |S] for the cold and hot runs respectively. 
For the cold run, G{f) calculated from the raw lightcurve is 
accurately approximated by a power law with an index of -2 



O 




Figure 6. Power spectra of the predicted emission from individ- 
ual annuli of the hot disc simulation. The system inclination is 
i = 50°. From bottom upward, we show results for annuli defined 
by 6 < r/rms < 7, 8 < r/rms < 9, and 10 < r/rms < 11. 
Below a break frequency, which is comparable to the orbital 
frequency, the power is flat or slowly declining (roughly brack- 
eted by G{f) = const and G(/) oc At high frequencies, 
G(/) oc !/-3'5. 



for all frequencies which are accessible using the simulation. 
There are no significant changes in the slope of G{f) with 
increasing system inclination. Figure |7| also shows power 
spectra computed from the NT-normalized lightcurves. Al- 
though not identical, these are almost indistinguishable from 
power spectra calculated using the raw lightcurve. We con- 
clude that the presence of torque at the marginally stable or- 
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0.1 1.0 0.1 1.0 

f f 

Figure 7. Upper panels: Power spectra G(/) for the raw (left) and NT- normalized (right) light curves resulting from the cold simulation. 
The frequency / has been rescaled such that the orbital frequency at the radius of marginal stability corresponds to / = 1. Shown here 
are the power spectra for i = 1°, i = 30° (offset by 10^), i = 50° (offset by 10^), and i = 80° (offset by 10^). Lower panels: The ratio 
G(/, i) /G{f, i = 1) for i = 50° and i = 80° . This representation helps isolate the effect of inclination on the power spectrum. See text 
for discussion. 




0.1 1.0 0.1 1.0 

f f 



Figure 8. Same as the upper panels in Fig.|7] except here computed for the hot simulation. 
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bit - at least at the rather modest levels obtained in our cal- 
culations ~ does not significantly affect the predicted power 
spectrum. 

The power spectra shown in Figure Q appear to be ac- 
curately described by single power laws over the range of 
frequencies accessible using the simulation. Observationally, 
high frequency QPOs in black hole systems are typically 
weak features. For example, the broad QPO at 300 Hz re- 
ported by Remillard et al. (1999) in GRO J 1655-40 had an 
amplitude of less that 1%. Substantially longer simulations 
would be needed before we could measure the theoretical 
power spectrum accurately enough to detect broad QPOs of 
such small amplitude. We can, however, place limits on the 
possible existence of stronger QPOs that could arise from 
the 'lighthouse' mechanism discussed by Abramowicz et al. 
(1991) and Karas (1999). Such QPOs, since they originate 
from the combination of hotspots with relativistic beam- 
ing, should become stronger at higher system inclinations. 
To investigate this possibility, the lower panels of Figure 
show the relative power for inclined systems as compared to 
face-on systems. For i = 50°, any QPOs at frequencies com- 
parable to the orbital frequency contribute at most ~ 20% 
of the broadband power. 

Figure |H| shows the corresponding power spectra ex- 
tracted from the hot simulation. The results are generally 
consistent, though the deviations from a power law form for 
G{f) are larger. At higher inclinations there are some in- 
dications of the presence of excess power at a frequency of 
around 0.2-0.3 times the orbital frequency at Vms- Excess 
power is also seen in the cold simulation at similar frequen- 
cies, though at a lower level. Longer simulations are required 
to investigate the reality of such weak, broad features. 

Since the emission from a particular annulus has a 
power-spectrum that breaks at the orbital frequency of 
that annulus, we expect that the disk-integrated power- 
spectrum should possess a high-frequency break correspond- 
ing to the innermost emitting annulus. The overall power- 
spectra should therefore steepen appreciably above / ~ 1. 
In fact, we do not see such a high-frequency break in Fig- 
ures and |H1 We have investigated the source of this addi- 
tional high frequency noise in our lightcurves, and find that 
it originates from emission from the outermost parts of the 
simulated disc. Our suspicion is that the spatial discretiza- 
tion of the disc into cells - which are physically largest in 
the outer disc - introduces a low level of high frequency 
fiicker into the lightcurve. This hampers our ability to see 
high-frequency cutoffs in the disk-integrated power-spectra. 



5 IRON LINE VARIABILITY 

The broad iron line is one of the best understood obser- 
vational probes of black hole accretion discs at the present 
time (e.g., see reviews by Fabian et al. 2000; Reynolds & 
Nowak 2003). Thus, it is useful to examine the predictions 
of our simulation for temporal variability of the iron line. 

The predicted temporal variability of the iron line obvi- 
ously depends on how that line is excited. Different authors 
have computed the iron line response to flares in the inner 
regions of the flow, using a variety of different geometries 
and kinematics for the source of exciting photons (Fabian et 
al. 1989; Stella 1990, Matt & Perola 1992; Reynolds et al. 




4 5 6 7 

Energy 



Figure 10. Iron line intensity as a function of energy and time 
for a 20 orbit segment of our simulation, assuming an inclination 
of i = 20°. 

1999; Ruszkowski 2000). In this work, we assume that the 
local iron line intensity coming from the surface of the disc 
is proportional to the local continuum intensity (and, hence, 
is proportional to the local vertically integrated azimuthal 
magnetic stress). The line intensity at image-plane position 
{x,y), time t and energy E is then given by 

hnc{x,y,t;E)<xg''{x,y)F{r,(l),Umit)S(E-gEo), (20) 

where Eo is the rest-frame energy of the emission line and 
F{r, (f), t) is taken to be either the raw of NT-normalized disc 
flux. Note that the K-correction is not relevant for comput- 
ing line fluxes and, in this expression, one power of g results 
from the transformation of the frequency-element. The to- 
tal simulated observed line is obtained from integrating this 
expression over the image plane. Figure |5] shows instanta- 
neous iron line profile resulting from the hot simulation for 
an inclination of i = 40°. The skewed and redshifted pro- 
file is caused by a combination of the line-of-sight Doppler 
effect, the transverse Dopper effect (i.e. special relativistic 
time dilation) and gravitational redshifts (see Fabian et al. 
2000 and Reynolds & Nowak 2003 for a detailed discussion 
of the physics determining these line profiles). 

Under the assumptions employed above, the simulated 
iron line at any particular instant is a weighted redshift map 
of the turbulent emissivity pattern across the disc. Thus, 
the temporal and azimuthal variability caused by the MHD 
turbulence will produce a level of inevitable and irreducible 
temporal variability in the iron line profile. This can be seen 
by comparing the instantaneous line profiles in Figure |3| 
From this figure, it can be seen that the iron line from the 
hot disc displays dramatic (factor of ~ 2) variability between 
independent time slices, whereas the line profile from the 
cold disc is more stable, showing deviations at only the ~ 30 
percent level. 

The nature of the iron line variability is demonstrated 
more clearly in Figure ITUI which shows the intensity of the 
simulated iron line emission as a function of energy and time 
for a representative interval during our simulations, assum- 
ing an inclination of i = 20°. This is a plausible inclination 
for many Seyfert-1 nuclei, the systems in which broad iron 
lines have been best studied to date. One can see a tremen- 
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(a) Cold run 



(b) Hot run 




Energy (keV) Energy (keV) 

Figure 9. Instantaneous simulated iron lino profiles for the cold (left) and hot (right) simulations. Note that the line flux is plotted on 
a logarithmic vertical scale. 



dous degree of structure in the iron line profile as a function 
of both energy and time. Of particular note are the distinc- 
tive arcs that correspond to the orbital motion of especially 
bright regions within the turbulent flow, visible in Figure |31 
These arcs essentially trace out (in the energy-time plane) 
the paths of test particle with orbits close to the radius of 
marginally stability. If they could be detected, they would 
provide a quantitative test of the flow extremely close to the 
black hole. High temporal resolution is required to observe 
such features, but the required capabilities should be read- 
ily within reach of the planned Constellation- X and XEUS 
missions. 

In fact, observational evidence for transient non- 
axisymmetric structure already exists. C/iandra-HETG ob- 
servations of the Seyfert-1 galaxy NGC 3516 have found sub- 
structure within the well known broad iron line profile of 
this object (Turner et al. 2002). The substructure has the 
form of narrow line-line features at rest-frame energies of 
5.6 keV, 6.2 keV and 7.1 keV. Plausible explanations exist 
for the latter two features; the 7.1 keV feature can be read- 
ily attributed to the fluorescent K/3 emission of cold iron, 
and the 6.2 keV feature might be the 'Compton shoulder' 
of the prominent 6.4 keV Ka fluorescent line of cold iron. 
However, the 5.6 keV feature cannot be attributed to any 
expected features. Turner et al. suggest that this feature is, 
in fact, iron Ka fluorescence with a rest- frame energy of 
6.4 keV from an anomalously bright annulus or spot on the 
accretion disc surface. The flux contained in the 5.6 keV fea- 
ture is approximately 5% of that contained in the overall 
broad iron line. 

Is it possible to relate the substructure observed in the 
NGC 3516 iron line with the turbulence induced substruc- 
ture seen in our simulations? While large amplitude sub- 
structure is present in virtually every instantaneous simu- 
lated line profile, it is much diminished once the line profiles 
are averaged over more than a few orbits. The Chandra- 
HETG observation of Turner et al. that yielded the clear 
detection of the iron fine substructure was 8 X 10* s in du- 
ration. Assuming a black hole mass of 2.3 X lO'^ Mq (esti- 
mated using the velocity dispersion from di Nella et al. 1995, 
in conjunction with the scaling relations of Gebhardt et al. 



2000 or Ferrarese & Merritt 2000)"^, this corresponds to ap- 
proximately 20 orbital periods of the accretion disc at the 
radius of marginal stability. Thus, to compare our simula- 
tions to these data, we averaged our simulated line profiles 
over 20 orbits of the inner disc. We assumed a disc inclina- 
tion of i = 40° (Wu & Han 2001). Given the useful length of 
our simulation, five independent 20-orbit line profiles could 
be formed. We found that 4 out of 5 of these profiles pos- 
sessed no strong deviations from the average. In one profile, 
however, there was a slight excess centered at an energy of 
5.5 keV, possessing 2-3% of the total flux. Therefore, while 
the true interpretation of the Chandra data for NGC 3516 
remains unclear, turbulence-induced substructure is a viable 
possibility even when one considers that the disc has likely 
undergone several orbits during the observation. 



6 SUMMARY 

In this paper, we have explored the predicted variability 
of accretion discs around Schwarzschild black holes. The 
temporal and spatial fluctuations in the magnetic stress 
which drives accretion have been calculated using three- 
dimensional global simulations in a pseudo-Newtonian po- 
tential, which ought to provide a reasonable first approxi- 
mation to the disc dynamics outside the marginally stable 
orbit. By making the further critical assumption that the 
radiation from the disc follows the local stress, we have cal- 
culated the predicted emission that would be seen by an 
observer using a fully relativistic ray-tracing code that ac- 
counts for the relativistic effects of beaming, gravitational 
redshift, and light bending. 

Our results for the temporal power spectrum of the pre- 
dicted emission are consistent with those of previous authors 
(Kawaguchi et al. 2000; Hawley & Krolik 2001, 2002). Fluc- 
tuations in the magnetic stress, when integrated across the 
disc, produce a high frequency power spectrum described 

^ Independent estimates of the black hole mass are comparable. 
For example, Onken et al. (2003) estimate Mbh = 1-7 x lO'' Mq 
from reverberation mapping. 
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by a power law with slope /^^. To a first approximation, 
observers viewing the system at a wide range of inclinations 
- from face-on to almost edge-on - are predicted to derive 
similar power spectra. In more detail, high inclinations are 
found to boost the relative power at frequencies a factor of 
a few below the orbital frequency at the marginally stable 
orbit. There is little change in the predicted power spectra 
due to the presence or absence of modest levels of mag- 
netic stress at the marginally stable orbit. No strong QPOs 
are seen in the simulations, though weak features could be 
present and would not be clearly observed given the limited 
statistics available from « 10^ orbit simulations. 

The broad iron line provides an alternate probe of con- 
ditions in the inner disc, and has the advantage of retain- 
ing velocity information. We have used the simulations to 
quantify the predicted iron line variability within a model 
in which the local iron line emissivity tracks the locally gen- 
erated continuum emission. This model represents the oppo- 
site extreme from reverberation models, in which the source 
of line exciting flux is often assumed to be widely spatially 
separated from the disc. We find that order unity fluctu- 
ations in the line flux at fixed energy occur for discs with 
h/r ~ 0.1, with the fiuctuation amplitude increasing rapidly 
with increasing disc thickness. Averaging over at least several 
orbits at the radii of interest is needed before meaningful av- 
erages can be extracted (see also Papaloizou & Nelson 2003; 
Winters, Balbus & Hawley 2003). Observationally, this im- 
plies that long integrations, extending over tens of orbital 
timescales, are required in order to reduce the fluctuation 
amplitude to percent levels and thereby recover the mean 
line profile. 
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